section ‹Missing Matrix Operations›

text ‹In this theory we provide an operation that can change a single
  row in a matrix efficiently, and all other rows in the matrix implementation
  will be reused.›

(* TODO: move this part into JNF-AFP-entry *)

theory Matrix_Change_Row
  imports 
    Jordan_Normal_Form.Matrix_IArray_Impl
    Polynomial_Interpolation.Missing_Unsorted
begin

definition change_row :: "nat ⇒ (nat ⇒ 'a ⇒ 'a) ⇒ 'a mat ⇒ 'a mat" where
  "change_row k f A = mat (dim_row A) (dim_col A) (λ (i,j). 
     if i = k then f j (A $$ (k,j)) else A $$ (i,j))"

lemma change_row_carrier[simp]: 
  "(change_row k f A ∈ carrier_mat nr nc) = (A ∈ carrier_mat nr nc)" 
  "dim_row (change_row k f A) = dim_row A" 
  "dim_col (change_row k f A) = dim_col A" 
  unfolding change_row_def carrier_mat_def by auto

lemma change_row_index[simp]: "A ∈ carrier_mat nr nc ⟹ i < nr ⟹ j < nc ⟹
  change_row k f A $$ (i,j) = (if i = k then f j (A $$ (k,j)) else A $$ (i,j))" 
  "i < dim_row A ⟹ j < dim_col A ⟹ change_row k f A $$ (i,j) = (if i = k then f j (A $$ (k,j)) else A $$ (i,j))" 
  unfolding change_row_def by auto

lift_definition change_row_impl :: "nat ⇒ (nat ⇒ 'a ⇒ 'a) ⇒ 'a mat_impl ⇒ 'a mat_impl" is
  "λ k f (nr,nc,A). let Ak = IArray.sub A k; Arows = IArray.list_of A;
     Ak' = IArray.IArray (map (λ (i,c). f i c) (zip [0 ..< nc] (IArray.list_of Ak)));
     A' = IArray.IArray (Arows [k := Ak'])
     in (nr,nc,A')" 
proof (auto, goal_cases)
  case (1 k f nc b row)
  show ?case 
  proof (cases b)
    case (IArray rows)
    with 1 have "row ∈ set rows ∨ k < length rows 
       ∧ row = IArray (map (λ (i,c). f i c) (zip [0 ..< nc] (IArray.list_of (rows ! k))))"
      by (cases "k < length rows", auto simp: set_list_update dest: in_set_takeD in_set_dropD)
    with 1 IArray show ?thesis by (cases, auto)
  qed
qed

lemma change_row_code[code]: "change_row k f (mat_impl A) = (if k < dim_row_impl A 
  then mat_impl (change_row_impl k f A) 
  else Code.abort (STR ''index out of bounds in change_row'') (λ _. change_row k f (mat_impl A)))"
  (is "?l = ?r")
proof (cases "k < dim_row_impl A")
  case True
  hence id: "?r = mat_impl (change_row_impl k f A)" by simp
  show ?thesis unfolding id unfolding change_row_def
  proof (rule eq_matI, goal_cases)
    case (1 i j)
    thus ?case using True
      by (transfer, auto simp: mk_mat_def)
  qed (transfer, auto)+
qed simp

end

section ‹Signed Modulo Operation›

theory Signed_Modulo
  imports 
    Berlekamp_Zassenhaus.Poly_Mod
    Sqrt_Babylonian.Sqrt_Babylonian_Auxiliary
begin

text ‹The upcoming definition of symmetric modulo 
  is different to the HOL-Library-Signed\_Division.smod, since
  here the modulus will be in range $\{-m/2,...,m/2\}$, 
  whereas there -1 symmod m = m - 1.

  The advantage of have range $\{-m/2,...,m/2\}$ is that small negative
  numbers are represented by small numbers.

  One limitation is that the symmetric modulo is only working properly,
  if the modulus is a positive number.›

definition sym_mod :: "int ⇒ int ⇒ int" (infixl "symmod" 70) where
  "sym_mod x y = poly_mod.inv_M y (x mod y)"

lemma sym_mod_code[code]: "sym_mod x y = (let m = x mod y
   in if m + m ≤ y then m else m - y)" 
  unfolding sym_mod_def poly_mod.inv_M_def Let_def ..

lemma sym_mod_zero[simp]: "n symmod 0 = n" "n > 0 ⟹ 0 symmod n = 0"
  unfolding sym_mod_def poly_mod.inv_M_def by auto

lemma sym_mod_range: "y > 0 ⟹ x symmod y ∈ {- ((y - 1) div 2) .. y div 2}"
  unfolding sym_mod_def poly_mod.inv_M_def using pos_mod_bound[of y x]
  by (cases "x mod y ≥ y", auto) 
    (smt (verit) Euclidean_Division.pos_mod_bound Euclidean_Division.pos_mod_sign half_nonnegative_int_iff)+

text ‹The range is optimal in the sense that exactly y elements can be represented.›
lemma card_sym_mod_range: "y > 0 ⟹ card {- ((y - 1) div 2) .. y div 2} = y" 
  by simp

lemma sym_mod_abs: "y > 0 ⟹ ¦x symmod y¦ < y"
  "y ≥ 1 ⟹ ¦x symmod y¦ ≤ y div 2"
  using sym_mod_range[of y x] by auto


lemma sym_mod_sym_mod[simp]: "x symmod y symmod y = x symmod (y :: int)" 
  unfolding sym_mod_def using poly_mod.M_def poly_mod.M_inv_M_id by auto

lemma sym_mod_diff_eq: "(a symmod c - b symmod c) symmod c = (a - b) symmod c" 
  unfolding sym_mod_def
  by (metis mod_diff_cong mod_mod_trivial poly_mod.M_def poly_mod.M_inv_M_id)

lemma sym_mod_sym_mod_cancel: "c dvd b ⟹ a symmod b symmod c = a symmod c" 
  using mod_mod_cancel[of c b] unfolding sym_mod_def
  by (metis poly_mod.M_def poly_mod.M_inv_M_id)

lemma sym_mod_diff_right_eq: "(a - b symmod c) symmod c = (a - b) symmod c" 
  using sym_mod_diff_eq by (metis sym_mod_sym_mod)

lemma sym_mod_mult_right_eq: "a * (b symmod c) symmod c = a * b symmod c" 
  unfolding sym_mod_def by (metis poly_mod.M_def poly_mod.M_inv_M_id mod_mult_right_eq)

lemma dvd_imp_sym_mod_0 [simp]:
  "b symmod a = 0" if "a > 0" "a dvd b"
  unfolding sym_mod_def poly_mod.inv_M_def using that by simp

lemma sym_mod_0_imp_dvd [dest!]:
  "b dvd a" if "a symmod b = 0"
  using that unfolding sym_mod_def poly_mod.inv_M_def
  by (smt (verit) Euclidean_Division.pos_mod_bound dvd_eq_mod_eq_0)

definition sym_div :: "int ⇒ int ⇒ int" (infixl "symdiv" 70) where
  "sym_div x y = (let d = x div y; m = x mod y in 
       if m + m ≤ y then d else d + 1)"

lemma of_int_mod_integer: "(of_int (x mod y) :: integer) = (of_int x :: integer) mod (of_int y)" 
  using integer_of_int_eq_of_int modulo_integer.abs_eq by presburger

lemma sym_div_code[code]: 
  "sym_div x y = (let yy = integer_of_int y in 
     (case divmod_integer (integer_of_int x) yy
     of (d, m) ⇒ if m + m ≤ yy then int_of_integer d else (int_of_integer (d + 1))))"
  unfolding sym_div_def Let_def divmod_integer_def split
  apply (rule if_cong, subst of_int_le_iff[symmetric], unfold of_int_add)
  by (subst (1 2) of_int_mod_integer, auto)

lemma sym_mod_sym_div: assumes y: "y > 0" shows "x symmod y = x - sym_div x y * y"
proof -
  let ?z = "x - y * (x div y)" 
  let ?u = "y * (x div y)" 
  have "x = y * (x div y) + x mod y" using y by simp
  hence id: "x mod y = ?z" by linarith
  have "x symmod y = poly_mod.inv_M y ?z" unfolding sym_mod_def id by auto
  also have "… = (if ?z + ?z ≤ y then ?z else ?z - y)" unfolding poly_mod.inv_M_def ..
  also have "… = x - (if (x mod y) + (x mod y) ≤ y then x div y else x div y + 1) * y" 
    by (simp add: algebra_simps id)
  also have "(if (x mod y) + (x mod y) ≤ y then x div y else x div y + 1) = sym_div x y" 
    unfolding sym_div_def Let_def ..
  finally show ?thesis .
qed
  
lemma dvd_sym_div_mult_right [simp]:
  "(a symdiv b) * b = a" if "b > 0" "b dvd a"
  using sym_mod_sym_div[of b a] that by simp

lemma dvd_sym_div_mult_left [simp]:
  "b * (a symdiv b) = a" if "b > 0" "b dvd a"
  using dvd_sym_div_mult_right[OF that] by (simp add: ac_simps)


end

section ‹Storjohann's Lemma 13›

text ‹This theory contains the result that one can always perform a mod-operation on
  the entries of the $d\mu$-matrix.›

theory Storjohann_Mod_Operation
  imports 
    LLL_Basis_Reduction.LLL_Certification
    Signed_Modulo
begin 

lemma map_vec_map_vec: "map_vec f (map_vec g v) = map_vec (f o g) v" 
  by (intro eq_vecI, auto)

context semiring_hom
begin

(* TODO: move *)
lemma mat_hom_add: assumes A: "A ∈ carrier_mat nr nc" and B: "B ∈ carrier_mat nr nc"
  shows "mat⇩h (A + B) = mat⇩h A + mat⇩h B"
  by (intro eq_matI, insert A B, auto simp: hom_add)
end

text ‹We now start to prove lemma 13 of Storjohann's paper.›
context
  fixes A I :: "'a :: field mat" and n :: nat
  assumes A: "A ∈ carrier_mat n n" 
  and det: "det A ≠ 0" 
  and I: "I = the (mat_inverse A)" 
begin
lemma inverse_via_det: "I * A = 1⇩m n" "A * I = 1⇩m n" "I ∈ carrier_mat n n" 
  "I = mat n n (λ (i,j). det (replace_col A (unit_vec n j) i) / det A)"
proof -
  from det_non_zero_imp_unit[OF A det] 
  have Unit: "A ∈ Units (ring_mat TYPE('a) n n)" .
  from mat_inverse(1)[OF A, of n] Unit I have "mat_inverse A = Some I" 
    by (cases "mat_inverse A", auto)
  from mat_inverse(2)[OF A this]
  show left: "I * A = 1⇩m n" and right: "A * I = 1⇩m n" and I: "I ∈ carrier_mat n n" 
    by blast+
  {
    fix i j
    assume i: "i < n" and j: "j < n" 
    from I i j have cI: "col I j $ i = I $$ (i,j)" by simp
    from j have uv: "unit_vec n j ∈ carrier_vec n" by auto
    from j I have col: "col I j ∈ carrier_vec n" by auto
    from col_mult2[OF A I j, unfolded right] j
    have "A *⇩v col I j = unit_vec n j" by simp
    from cramer_lemma_mat[OF A col i, unfolded this cI]
    have "I $$ (i,j) = det (replace_col A (unit_vec n j) i) / det A" using det by simp
  }
  thus "I = mat n n (λ (i,j). det (replace_col A (unit_vec n j) i) / det A)"
    by (intro eq_matI, use I in auto)
qed

lemma matrix_for_singleton_entry: assumes i: "i < n" and 
  j: "j < n" 
  and Rdef: "R = mat n n ( λ ij. if ij = (i,j) then c :: 'a else 0)" 
shows "mat n n
   (λ(i', j'). if i' = i then c * det (replace_col A (unit_vec n j') j) / det A
       else 0) * A = R" 
proof -
  note I = inverse_via_det(3)
  have R: "R ∈ carrier_mat n n" unfolding Rdef by auto
  have "(R * I) * A = R * (I * A)" using I A R by auto
  also have "I * A = 1⇩m n" unfolding inverse_via_det(1) ..
  also have "R * … = R" using R by simp
  also have "R * I = mat n n (λ (i',j'). row R i' ∙ col I j')"
    using I R unfolding times_mat_def by simp
  also have "… = mat n n ( λ (i',j'). if i' = i then c * I $$ (j, j') else 0)" 
    (is "mat n n ?f = mat n n ?g")
  proof -
    {
      fix i' j'
      assume i': "i' < n" and j': "j' < n" 
      have "?f (i',j') = ?g (i',j')" 
      proof (cases "i' = i")
        case False
        hence "row R i' = 0⇩v n" unfolding Rdef using i'
          by (intro eq_vecI, auto simp: Matrix.row_def)
        thus ?thesis using False i' j' I by simp
      next
        case True
        hence "row R i' = c ⋅⇩v unit_vec n j" unfolding Rdef using i' j' i j
          by (intro eq_vecI, auto simp: Matrix.row_def)
        with True show ?thesis using i' j' I j by simp
      qed
    }
    thus ?thesis by auto
  qed
  finally show ?thesis unfolding inverse_via_det(4) using j 
    by (auto intro!: arg_cong[of _ _ "λ x. x * A"])
qed
end

lemma (in gram_schmidt_fs_Rn) det_M_1: "det (M m) = 1" 
proof -
  have "det (M m) = prod_list (diag_mat (M m))" 
    by (rule det_lower_triangular[of m], auto simp: μ.simps)
  also have "… = 1" 
    by (rule prod_list_neutral, auto simp: diag_mat_def μ.simps)
  finally show ?thesis .
qed

context gram_schmidt_fs_int
begin
lemma assumes IM: "IM = the (mat_inverse (M m))" 
  shows inv_mu_lower_triangular: "⋀ k i. k < i ⟹ i < m ⟹ IM $$ (k, i) = 0"
  and inv_mu_diag: "⋀ k. k < m ⟹ IM $$ (k, k) = 1"
  and d_inv_mu_integer: "⋀ i j. i < m ⟹ j < m ⟹ d i * IM $$ (i,j) ∈ ℤ" 
  and inv_mu_inverse: "IM * M m = 1⇩m m" "M m * IM = 1⇩m m" "IM ∈ carrier_mat m m" 
proof - 
  note * = inverse_via_det[OF M_dim(3) _ IM, unfolded det_M_1]
  from * show inv: "IM * M m = 1⇩m m" "M m * IM = 1⇩m m" 
    and IM: "IM ∈ carrier_mat m m"  by auto
  from * have IM_det: "IM = mat m m (λ(i, j). det (replace_col (M m) ((unit_vec m) j) i))" 
    by auto
  from matrix_equality have "IM * FF = IM * ((M m) * Fs)" by simp
  also have "… = (IM * M m) * Fs" using M_dim(3) IM Fs_dim(3)
    by (metis assoc_mult_mat)
  also have "… = Fs" unfolding inv using Fs_dim(3) by simp
  finally have equality: "IM * FF = Fs" .
  {
    fix i k
    assume i: "k < i" "i < m" 
    show "IM $$ (k, i) = 0" using i M_dim unfolding IM_det
      by (simp, subst det_lower_triangular[of m], auto simp: replace_col_def μ.simps diag_mat_def)
  } note IM_lower_triag = this
  {
    fix k
    assume k: "k < m" 
    show "IM $$ (k,k) = 1" using k M_dim unfolding IM_det
      by (simp, subst det_lower_triangular[of m], auto simp: replace_col_def μ.simps diag_mat_def
        intro!: prod_list_neutral)
  } note IM_diag_1 = this
  {
    fix k
    assume k: "k < m" 
    let ?f = "λ i. IM $$ (k, i) ⋅⇩v fs ! i" 
    let ?sum = "M.sumlist (map ?f [0..<m])" 
    let ?sumk = "M.sumlist (map ?f [0..<k])" 
    have set: "set (map ?f [0..<m]) ⊆ carrier_vec n" using fs_carrier by auto
    hence sum: "?sum ∈ carrier_vec n" by simp
    from set k have setk: "set (map ?f [0..<k]) ⊆ carrier_vec n" by auto
    hence sumk: "?sumk ∈ carrier_vec n" by simp
    from sum have dim_sum: "dim_vec ?sum = n" by simp
    have "gso k = row Fs k" using k by auto
    also have "… = row (IM * FF) k" unfolding equality ..
    also have "IM * FF = mat m n (λ (i,j). row IM i  ∙ col FF j)" 
      unfolding times_mat_def using IM FF_dim by auto
    also have "row … k = vec n (λ j. row IM k ∙ col FF j)" 
      unfolding Matrix.row_def using IM FF_dim k by auto
    also have "… = vec n (λ j. ∑ i < m. IM $$ (k, i) * fs ! i $ j)" 
      by (intro eq_vecI, insert IM k, auto simp: scalar_prod_def Matrix.row_def intro!: sum.cong)
    also have "… = ?sum" 
      by (intro eq_vecI, insert IM, unfold dim_sum, subst sumlist_vec_index, 
        auto simp: o_def sum_list_sum_nth intro!: sum.cong)
    also have "[0..<m] = [0..<k] @ [k] @ [Suc k ..<m]" using k
      by (simp add: list_trisect)
    also have "M.sumlist (map ?f …) = ?sumk + 
      (?f k + M.sumlist (map ?f [Suc k ..< m]))" 
      unfolding map_append 
      by (subst M.sumlist_append; (subst M.sumlist_append)?, insert k fs_carrier, auto)
    also have "M.sumlist (map ?f [Suc k ..< m]) = 0⇩v n" 
      by (rule sumlist_neutral, insert IM_lower_triag, auto)
    also have "IM $$ (k,k) = 1" using IM_diag_1[OF k] .
    finally have gso: "gso k = ?sumk + fs ! k"  using k by simp
    define b where "b = vec k (λ j. fs ! j ∙ fs ! k)" 
    {
      fix j
      assume jk: "j < k" 
      with k have j: "j < m" by auto
      have "fs ! j ∙ gso k = fs ! j ∙ (?sumk + fs ! k)" 
        unfolding gso by simp
      also have "fs ! j ∙ gso k = 0" using jk k
        by (simp add: fi_scalar_prod_gso gram_schmidt_fs.μ.simps)
      also have "fs ! j ∙ (?sumk + fs ! k)
         = fs ! j ∙ ?sumk + fs ! j ∙ fs ! k" 
        by (rule scalar_prod_add_distrib[OF _ sumk], insert j k, auto)
      also have "fs ! j ∙ fs ! k = b $ j" unfolding b_def using jk by simp
      finally have "b $ j = - (fs ! j ∙ ?sumk)" by linarith
    } note b_index = this
    let ?x = "vec k (λ i. - IM $$ (k, i))" 
    have x: "?x ∈ carrier_vec k" by auto
    from k have km: "k ≤ m" by simp 
    have bGx: "b = Gramian_matrix fs k *⇩v (vec k (λ i. - IM $$ (k, i)))" 
      unfolding Gramian_matrix_alt_alt_def[OF km]
    proof (rule eq_vecI; simp)
      fix i
      assume i: "i < k" 
      have "b $ i = - (∑x←[0..<k]. fs ! i ∙ (IM $$ (k, x) ⋅⇩v fs ! x))" 
        unfolding b_index[OF i]
        by (subst scalar_prod_right_sum_distrib, insert setk i k, auto simp: o_def)
      also have "… = vec k (λj. fs ! i ∙ fs ! j) ∙ vec k (λi. - IM $$ (k, i))" 
        by (subst (3) scalar_prod_def, insert i k, auto simp: o_def sum_list_sum_nth simp flip: sum_negf
          intro!: sum.cong)
      finally show "b $ i = vec k (λj. fs ! i ∙ fs ! j) ∙ vec k (λi. - IM $$ (k, i))" .
    qed (simp add: b_def)
    have G: "Gramian_matrix fs k ∈ carrier_mat k k" 
      unfolding Gramian_matrix_alt_alt_def[OF km] by simp
    from cramer_lemma_mat[OF G x, folded bGx Gramian_determinant_def]
    have "i < k ⟹ 
      d k * IM $$ (k, i) = - det (replace_col (Gramian_matrix fs k) (vec k (λ j. fs ! j ∙ fs ! k)) i)" 
      for i unfolding b_def by simp
  } note IM_lower_values = this
  {
    fix i j
    assume i: "i < m" and j: "j < m" 
    from i have im: "i ≤ m" by auto
    consider (1) "j < i" | (2) "j = i" | (3) "i < j" by linarith
    thus "d i * IM $$ (i,j) ∈ ℤ"
    proof cases
      case 1
      show ?thesis unfolding IM_lower_values[OF i 1] replace_col_def Gramian_matrix_alt_alt_def[OF im]
        by (intro Ints_minus Ints_det, insert i j, auto intro!: Ints_scalar_prod[of _ n] fs_int)
    next
      case 3
      show ?thesis unfolding IM_lower_triag[OF 3 j] by simp
    next
      case 2
      show ?thesis unfolding IM_diag_1[OF i] 2 using i unfolding Gramian_determinant_def
         Gramian_matrix_alt_alt_def[OF im]
        by (intro Ints_mult Ints_det, insert i j, auto intro!: Ints_scalar_prod[of _ n] fs_int)
    qed 
  }
qed

definition inv_mu_ij_mat :: "nat ⇒ nat ⇒ int ⇒ int mat" where
 "inv_mu_ij_mat i j c = (let
    B = mat m m (λ ij. if ij = (i,j) then c else 0);
    C = mat m m (λ (i,j). the_inv (of_int :: _ ⇒ 'a) (d i * the (mat_inverse (M m)) $$ (i,j)))
   in B * C + 1⇩m m)" 

lemma inv_mu_ij_mat: assumes i: "i < m" and ji: "j < i" 
  shows 
(* Effect on μ *)
   "map_mat of_int (inv_mu_ij_mat i j c) * M m =
    mat m m (λij. if ij = (i, j) then of_int c * d j else 0) + M m" (* only change value of μ_ij *)
(* Effect on A *)
  "A ∈ carrier_mat m n ⟹ c mod p = 0 ⟹ map_mat (λ x. x mod p) (inv_mu_ij_mat i j c * A) = 
    (map_mat (λ x. x mod p) A)" (* no change (mod p) *)
(* The transformation-matrix is ... *)
  "inv_mu_ij_mat i j c ∈ carrier_mat m m" (* ... of dimension m*m *)
  "i' < j' ⟹ j' < m ⟹ inv_mu_ij_mat i j c $$ (i',j') = 0" (* ... lower triangular *)
  "k < m ⟹ inv_mu_ij_mat i j c $$ (k,k) = 1" (* ... with diagonal all 1 *)  
proof -
  obtain IM where IM: "IM = the (mat_inverse (M m))" by auto
  let ?oi = "of_int :: _ ⇒ 'a" 
  let ?C = "mat m m (λ ij. if ij = (i,j) then ?oi c else 0)" 
  let ?D = "mat m m (λ (i,j). d i * IM $$ (i,j))" 
  have oi: "inj ?oi" unfolding inj_on_def by auto
  have C: "?C ∈ carrier_mat m m" by auto
  from i ji have j: "j < m" by auto
  from j have jm: "{0..<m} = {0..<j} ∪ {j} ∪ {Suc j..<m}" by auto
  note IM_props = d_inv_mu_integer[OF IM] inv_mu_inverse[OF IM]
  have mat_oi: "map_mat ?oi (inv_mu_ij_mat i j c) = ?C * ?D + 1⇩m m" (is "?MM = _")
    unfolding inv_mu_ij_mat_def Let_def IM[symmetric]
    apply (subst of_int_hom.mat_hom_add, force, force)
    apply (rule arg_cong2[of _ _ _ _ "(+)"])
     apply (subst of_int_hom.mat_hom_mult, force, force)
     apply (rule arg_cong2[of _ _ _ _ "(*)"])
      apply force
     apply (rule eq_matI, (auto)[3], goal_cases)
  proof -
    case (1 i j)
    from IM_props(1)[OF 1]
    show ?case unfolding Ints_def using the_inv_f_f[OF oi] by auto
  qed auto
  have "map_mat ?oi (inv_mu_ij_mat i j c) * M m = (?C * ?D) * M m + M m" unfolding mat_oi
    by (subst add_mult_distrib_mat[of _ m m], auto)
  also have "(?C * ?D) * M m = ?C * (?D * M m)" 
    by (rule assoc_mult_mat, auto)
  also have "?D = mat m m (λ (i,j). if i = j then d j else 0) * IM" (is "_ = ?E * _")
  proof (rule eq_matI, insert IM_props(4), auto simp: scalar_prod_def, goal_cases)
    case (1 i j)
    hence id: "{0..<m} = {0..<i} ∪ {i} ∪ {Suc i ..<m}" 
      by (auto simp add: list_trisect)
    show ?case unfolding id
      by (auto simp: sum.union_disjoint)
  qed
  also have "… * M m = ?E * (IM * M m)" 
    by (rule assoc_mult_mat[of _ m m], insert IM_props, auto)
  also have "IM * M m = 1⇩m m" by fact
  also have "?E * 1⇩m m = ?E" by simp
  also have "?C * ?E = mat m m (λ ij. if ij = (i,j) then ?oi c * d j else 0)" 
    by (rule eq_matI, auto simp: scalar_prod_def, auto simp: jm sum.union_disjoint)
  finally show "map_mat ?oi (inv_mu_ij_mat i j c) * M m = 
    mat m m (λ ij. if ij = (i,j) then ?oi c * d j else 0) + M m" .
  show carr: "inv_mu_ij_mat i j c ∈ carrier_mat m m"
    unfolding inv_mu_ij_mat_def by auto
  {
    assume k: "k < m" 
    have "of_int (inv_mu_ij_mat i j c $$ (k,k)) = ?MM $$ (k,k)" 
      using carr k by auto
    also have "… = (?C * ?D) $$ (k,k) + 1" unfolding mat_oi using k by simp
    also have "(?C * ?D) $$ (k,k) = 0" using k
      by (auto simp: scalar_prod_def, auto simp: jm sum.union_disjoint 
        inv_mu_lower_triangular[OF IM ji i])
    finally show "inv_mu_ij_mat i j c $$ (k,k) = 1" by simp
  }
  {
    assume ij': "i' < j'" "j' < m"  
    have "of_int (inv_mu_ij_mat i j c $$ (i',j')) = ?MM $$ (i',j')" 
      using carr ij' by auto
    also have "… = (?C * ?D) $$ (i',j')" unfolding mat_oi using ij' by simp
    also have "(?C * ?D) $$ (i',j') = (if i' = i then ?oi c * (d j * IM $$ (j, j')) else 0)" 
      using ij' i j by (auto simp: scalar_prod_def, auto simp: jm sum.union_disjoint)
    also have "… = 0" using inv_mu_lower_triangular[OF IM _ ij'(2), of j] ij' i ji by auto
    finally show "inv_mu_ij_mat i j c $$ (i',j') = 0" by simp
  }
  {
    assume A: "A ∈ carrier_mat m n" and c: "c mod p = 0" 
    let ?mod = "map_mat (λ x. x mod p)" 
    let ?C = "mat m m (λ ij. if ij = (i,j) then c else 0)" 
    let ?D = "mat m m (λ ij. if ij = (i,j) then 1 else (0 :: int))" 
    define B where "B = mat m m (λ (i,j). the_inv ?oi (d i * the (mat_inverse (M m)) $$ (i,j)))" 
    have B: "B ∈ carrier_mat m m" unfolding B_def by auto
    define BA where "BA = B * A" 
    have BA: "BA ∈ carrier_mat m n" unfolding BA_def using A B by auto
    define DBA where "DBA = ?D * BA" 
    have DBA: "DBA ∈ carrier_mat m n" unfolding DBA_def using BA by auto
    have "?mod (inv_mu_ij_mat i j c * A) = 
     ?mod ((?C * B + 1⇩m m) * A)" 
      unfolding inv_mu_ij_mat_def B_def by simp
    also have "(?C * B + 1⇩m m) * A = ?C * B * A + A" 
      by (subst add_mult_distrib_mat, insert A B, auto)
    also have "?C * B * A = ?C * BA" 
      unfolding BA_def
      by (rule assoc_mult_mat, insert A B, auto)
    also have "?C = c ⋅⇩m ?D" 
      by (rule eq_matI, auto)
    also have "… * BA = c ⋅⇩m DBA" using BA unfolding DBA_def by auto
    also have "?mod (… + A) = ?mod A" 
      by (rule eq_matI, insert DBA A c, auto simp: mult.assoc) 
    finally show "?mod (inv_mu_ij_mat i j c * A) = ?mod A" .
  }
qed   
end
 
lemma Gramian_determinant_of_int: assumes fs: "set fs ⊆ carrier_vec n" 
  and j: "j ≤ length fs" 
shows "of_int (gram_schmidt.Gramian_determinant n fs j)
  = gram_schmidt.Gramian_determinant n (map (map_vec rat_of_int) fs) j" 
proof -
  from j have j: "k < j ⟹ k < length fs" for k by auto
  show ?thesis
  unfolding gram_schmidt.Gramian_determinant_def
  by (subst of_int_hom.hom_det[symmetric], rule arg_cong[of _ _ det],
      unfold gram_schmidt.Gramian_matrix_def Let_def, subst of_int_hom.mat_hom_mult, force, force,
      unfold map_mat_transpose[symmetric],
      rule arg_cong2[of _ _ _ _ "λ x y. x * y⇧T"], insert fs[unfolded set_conv_nth] 
      j, (fastforce intro!: eq_matI)+)
qed

context LLL
begin

(* this lemma might also be useful for swap/add-operation *)
lemma multiply_invertible_mat: assumes lin: "lin_indep fs" 
  and len: "length fs = m" 
  and A: "A ∈ carrier_mat m m" 
  and A_invertible: "∃ B. B ∈ carrier_mat m m ∧ B * A = 1⇩m m" 
  and fs'_prod: "fs' = Matrix.rows (A * mat_of_rows n fs)" 
shows "lattice_of fs' = lattice_of fs" 
  "lin_indep fs'" 
  "length fs' = m" 
proof -
  let ?Mfs = "mat_of_rows n fs" 
  let ?Mfs' = "mat_of_rows n fs'" 
  from A_invertible obtain B where B: "B ∈ carrier_mat m m" and inv: "B * A = 1⇩m m" by auto
  from lin have fs: "set fs ⊆ carrier_vec n" unfolding gs.lin_indpt_list_def by auto
  with len have Mfs: "?Mfs ∈ carrier_mat m n" by auto
  from A Mfs have prod: "A * ?Mfs ∈ carrier_mat m n" by auto
  hence fs': "length fs' = m" "set fs' ⊆ carrier_vec n" unfolding fs'_prod
    by (auto simp: Matrix.rows_def Matrix.row_def)  
  have Mfs_prod': "?Mfs' = A * ?Mfs" 
    unfolding arg_cong[OF fs'_prod, of "mat_of_rows n"]
    by (intro eq_matI, auto simp: mat_of_rows_def)
  have "B * ?Mfs' = B * (A * ?Mfs)" 
    unfolding Mfs_prod' by simp
  also have "… = (B * A) * ?Mfs"
    by (subst assoc_mult_mat[OF _ A Mfs], insert B, auto)
  also have "B * A = 1⇩m m" by fact
  also have "… * ?Mfs = ?Mfs" using Mfs by auto
  finally have Mfs_prod: "?Mfs = B * ?Mfs'" ..  
  interpret LLL: LLL_with_assms n m fs 2
    by (unfold_locales, auto simp: len lin)
  from LLL.LLL_change_basis[OF fs'(2,1) B A Mfs_prod Mfs_prod']
  show latt': "lattice_of fs' = lattice_of fs" and lin': "gs.lin_indpt_list (RAT fs')" 
    and len': "length fs' = m" 
    by (auto simp add: LLL_with_assms_def)
qed

text ‹This is the key lemma.›
lemma change_single_element: assumes lin: "lin_indep fs" 
  and len: "length fs = m" 
  and i: "i < m" and ji: "j < i"  
  and A: "A = gram_schmidt_fs_int.inv_mu_ij_mat n (RAT fs)"    ― ‹the transformation matrix A›
  and fs'_prod: "fs' = Matrix.rows (A i j c * mat_of_rows n fs)" ― ‹fs' is the new basis›
  and latt: "lattice_of fs = L" 
shows "lattice_of fs' = L"
  "c mod p = 0 ⟹ map (map_vec (λ x. x mod p)) fs' = map (map_vec (λ x. x mod p)) fs" 
  "lin_indep fs'" 
  "length fs' = m" 
  "⋀ k. k < m ⟹ gso fs' k = gso fs k" 
  "⋀ k. k ≤ m ⟹ d fs' k = d fs k" 
  "i' < m ⟹ j' < m ⟹ 
    μ fs' i' j' = (if (i',j') = (i,j) then rat_of_int (c * d fs j) + μ fs i' j' else μ fs i' j')" 
  "i' < m ⟹ j' < m ⟹ 
    dμ fs' i' j' = (if (i',j') = (i,j) then c * d fs j * d fs (Suc j) + dμ fs i' j' else dμ fs i' j')" 
proof -
  let ?A = "A i j c" 
  let ?Mfs = "mat_of_rows n fs" 
  let ?Mfs' = "mat_of_rows n fs'" 
  from lin have fs: "set fs ⊆ carrier_vec n" unfolding gs.lin_indpt_list_def by auto
  with len have Mfs: "?Mfs ∈ carrier_mat m n" by auto
  interpret gsi: gram_schmidt_fs_int n "RAT fs"
    rewrites "gsi.inv_mu_ij_mat = A" using lin unfolding A
    by (unfold_locales, insert lin[unfolded gs.lin_indpt_list_def], auto simp: set_conv_nth)
  note A = gsi.inv_mu_ij_mat[unfolded length_map len, OF i ji, where c = c]
  from A(3) Mfs have prod: "?A * ?Mfs ∈ carrier_mat m n" by auto
  hence fs': "length fs' = m" "set fs' ⊆ carrier_vec n" unfolding fs'_prod
    by (auto simp: Matrix.rows_def Matrix.row_def)  
  have Mfs_prod': "?Mfs' = ?A * ?Mfs" 
    unfolding arg_cong[OF fs'_prod, of "mat_of_rows n"]
    by (intro eq_matI, auto simp: mat_of_rows_def)
  have detA: "det ?A = 1" 
    by (subst det_lower_triangular[OF A(4) A(3)], insert A, auto intro!: prod_list_neutral 
      simp: diag_mat_def)
  have "∃ B. B ∈ carrier_mat m m ∧ B * ?A = 1⇩m m" 
    by (intro exI[of _ "adj_mat ?A"], insert adj_mat[OF A(3)], auto simp: detA)
  from multiply_invertible_mat[OF lin len A(3) this fs'_prod] latt
  show latt': "lattice_of fs' = L" and lin': "gs.lin_indpt_list (RAT fs')" 
    and len': "length fs' = m" by auto
  interpret LLL: LLL_with_assms n m fs 2
    by (unfold_locales, auto simp: len lin)
  interpret fs: fs_int_indpt n fs
    by (standard, auto simp: lin)
  interpret fs': fs_int_indpt n fs'
    by (standard, auto simp: lin')
  {
    assume c: "c mod p = 0" 
    have id: "rows (map_mat f A) = map (map_vec f) (rows A)" for f A
      unfolding rows_def by auto
    have rows_id: "set fs ⊆ carrier_vec n ⟹ rows (mat_of_rows n fs) = fs" for fs
      unfolding mat_of_rows_def rows_def 
      by (force simp: Matrix.row_def set_conv_nth intro!: nth_equalityI)
    from A(2)[OF Mfs c]
    have "rows (map_mat (λx. x mod p) ?Mfs') = rows (map_mat (λx. x mod p) ?Mfs)" unfolding Mfs_prod'
      by simp
    from this[unfolded id rows_id[OF fs] rows_id[OF fs'(2)]]
    show "map (map_vec (λ x. x mod p)) fs' = map (map_vec (λ x. x mod p)) fs" .
  }
  {
    define B where "B = ?A" 
    have gs_eq: "k < m ⟹ gso fs' k = gso fs k" for k
    proof(induct rule: nat_less_induct)
      case (1 k)
      then show ?case 
      proof(cases "k = 0")
        case True
        then show ?thesis 
        proof -
          have "row ?Mfs' 0 = row ?Mfs 0"
          proof -
            have 2: "0∈ {0..<m}" and 3: "{1..<m} = {0..<m} - {0}" 
              and 4: "finite {0..<m}" using 1 by auto
            have "row ?Mfs' 0 = vec n (λj. row B 0 ∙ col ?Mfs j)" 
              using row_mult A(3) Mfs 1 Mfs_prod' unfolding B_def by simp
            also have "… = vec n (λj. (∑l∈{0..<m}. B $$ (0, l) * ?Mfs $$ (l, j)))"
              using Mfs A(3) len 1 B_def unfolding scalar_prod_def by auto
            also have "… = vec n (λj. B $$ (0, 0) * ?Mfs $$ (0, j) + 
              (∑l∈{1..<m}. B $$ (0, l) * ?Mfs $$ (l, j)))"
              using Groups_Big.comm_monoid_add_class.sum.remove[OF 4 2] 3
              by (simp add: ‹⋀g. sum g {0..<m} = g 0 + sum g ({0..<m} - {0})›)
            also have "… = row ?Mfs 0" 
              using A(4-) 1 unfolding B_def[symmetric] by (simp add: row_def)
            finally show ?thesis by (simp add: B_def Mfs_prod')
          qed
          then show ?thesis using True 1 fs'.f_carrier fs.f_carrier 
            fs'.gs.fs0_gso0 len' len gsi.fs0_gso0 by auto
        qed
      next
        case False
        then show ?thesis
        proof -
          have gso0kcarr: "gsi.gso ` {0 ..<k} ⊆ carrier_vec n"
            using 1(2) gsi.gso_carrier len by auto
          hence gsospancarr: "gs.span(gsi.gso ` {0 ..<k}) ⊆ carrier_vec n " 
            using span_is_subset2 by auto

          have fs'_gs_diff_span: 
            "(RAT fs') !  k - fs'.gs.gso k ∈ gs.span (gsi.gso ` {0 ..<k})"
          proof -
            define gs'sum where "gs'sum =
              gs.M.sumlist (map (λja. fs'.gs.μ k ja ⋅⇩v fs'.gs.gso ja) [0..<k])"
            define gssum where "gssum = 
              gs.M.sumlist (map (λja. fs'.gs.μ k ja ⋅⇩v gsi.gso ja) [0..<k])"
            have "set (map (λja. fs'.gs.μ k ja ⋅⇩v gsi.gso ja) [0..<k]) 
              ⊆ gs.span(gsi.gso ` {0 ..<k})" using 1(2) gs.span_mem gso0kcarr
              by auto
            hence gssumspan: "gssum ∈ gs.span(gsi.gso ` {0 ..<k})"
              using atLeastLessThan_iff gso0kcarr imageE set_map set_upt 
                vec_space.sumlist_in_span 
              unfolding gssum_def by (smt subsetD)
            hence gssumcarr: "gssum ∈ carrier_vec n" 
              using gsospancarr gssum_def by blast
            have sumid: "gs'sum = gssum"
            proof -
              have "map (λja. fs'.gs.μ k ja ⋅⇩v fs'.gs.gso ja) [0..<k] =
                map (λja. fs'.gs.μ k ja ⋅⇩v gsi.gso ja) [0..<k]"
                using 1 by simp
              thus ?thesis unfolding gs'sum_def gssum_def by argo
            qed
            have "(RAT fs') !  k = fs'.gs.gso k + gssum" 
              using fs'.gs.fs_by_gso_def len' False 1 sumid 
              unfolding gs'sum_def by auto
            hence "(RAT fs') !  k - fs'.gs.gso k = gssum" 
              using gssumcarr 1(2) len' by auto
            thus ?thesis using gssumspan by simp
          qed

          define v2 where "v2 = sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..< k])"
          have v2carr: "v2 ∈ carrier_vec n" 
          proof -
            have "set (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..< k]) ⊆ carrier_vec n"
              using len 1(2) fs.f_carrier by auto
            thus ?thesis unfolding v2_def by simp
          qed
          define ratv2 where "ratv2 = (map_vec rat_of_int v2)"
          have ratv2carr: "ratv2 ∈ carrier_vec n" 
            unfolding ratv2_def using v2carr by simp
          have fs'id: "(RAT fs') ! k = (RAT fs) ! k + ratv2"
          proof -
            have zkm: "[0..<m] = [0..<(Suc k)] @ [(Suc k)..<m]" using 1(2) 
              by (metis Suc_lessI append_Nil2 upt_append upt_rec zero_less_Suc)
            have prep: "set (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..<m]) ⊆ carrier_vec n" 
              using len fs.f_carrier by auto

            have "fs' ! k = vec n (λj. row B k ∙ col ?Mfs j)"
              using 1(2) Mfs B_def A(3) fs'_prod by simp
            also have "… = sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..<m])"
            proof -
              {
                fix i
                assume i: "i < n"
                have "(vec n (λj. row B k ∙ col ?Mfs j)) $ i = row B k ∙ col ?Mfs i" 
                  using i by auto
                also have "… = (∑j = 0..<m. B $$ (k, j) * ?Mfs $$ (j,i))" 
                  using A(3) unfolding B_def[symmetric] 
                  by (smt 1(2) Mfs R.finsum_cong' i atLeastLessThan_iff carrier_matD
                      dim_col index_col index_row(1) scalar_prod_def)
                also have "… = (∑j = 0..<m. B $$ (k, j) * (fs ! j $ i))"
                  by (metis (no_types, lifting) R.finsum_cong' atLeastLessThan_iff i
                      len mat_of_rows_index)
                also have "… = 
                  (∑j = 0..<m. (map (λja.  B $$ (k, ja) ⋅⇩v fs ! ja) [0..<m]) ! j $ i)"
                proof -
                  have "∀j<m. ∀i<n. B $$ (k, j) * (fs ! j $ i) = 
                    (map (λja.  B $$ (k, ja) ⋅⇩v fs ! ja) [0..<m]) ! j $ i" 
                    using 1(2) i A(3) len fs.f_carrier
                    unfolding B_def[symmetric] by auto
                  then show ?thesis using i by auto
                qed
                also have "… = sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..<m]) $ i"
                  using sumlist_nth i fs.f_carrier carrier_vecD len by simp
                finally have "(vec n (λj. row B k ∙ col ?Mfs j)) $ i =
                  sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..<m]) $ i" by auto
              }
              then show ?thesis using fs.f_carrier len dim_sumlist by auto
            qed
            also have "… = sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) 
              ([0..<(Suc k)] @ [(Suc k)..<m]))" 
              using zkm by simp
            also have "… = sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..<(Suc k)]) +
              sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [(Suc k)..<m])"
              (is "… = ?L2 + ?L3")
              using fs.f_carrier len dim_sumlist sumlist_append prep zkm by auto
            also have "?L3 = 0⇩v n"
              using A(4) fs.f_carrier len sumlist_nth carrier_vecD sumlist_carrier 
                prep zkm unfolding B_def[symmetric] by auto
            also have "?L2 = sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..<k]) +
              B $$ (k, k) ⋅⇩v fs ! k" using prep zkm sumlist_snoc by simp
            also have "… = sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..<k]) + fs ! k"
              using A(5) 1(2) unfolding B_def[symmetric] by simp
            finally have "fs' ! k = fs ! k + 
              sumlist (map (λja. B $$ (k, ja) ⋅⇩v fs ! ja) [0..<k])"
              using prep zkm by (simp add: M.add.m_comm)
            then have "fs' !  k = fs !  k + v2" unfolding v2_def by simp
            then show ?thesis using v2carr 1(2) len len' ratv2_def by force
          qed
          have ratv2span: "ratv2 ∈ gs.span (gsi.gso ` {0 ..<k})" 
          proof -
            have rat: "ratv2 = gs.M.sumlist
              (map (λj. of_int (B $$ (k, j)) ⋅⇩v (RAT fs) ! j) [0..<k])"
            proof -
              have "set (map (λj. of_int (B $$ (k, j)) ⋅⇩v (RAT fs) ! j) [0..<k]) 
                ⊆ carrier_vec n"
                using fs.f_carrier 1(2) len by auto
              hence carr: "gs.M.sumlist 
                (map (λj. of_int (B $$ (k, j)) ⋅⇩v (RAT fs) ! j) [0..<k]) ∈ carrier_vec n"
                by auto
              have "set (map (λj. B $$ (k, j) ⋅⇩v fs ! j) [0..<k]) ⊆ carrier_vec n"
                using fs.f_carrier 1(2) len by auto
              hence "⋀i j. i < n ⟹ j < k ⟹ of_int ((B $$ (k, j) ⋅⇩v fs ! j) $ i)
                = (of_int (B $$ (k, j)) ⋅⇩v (RAT fs) ! j) $ i"
                using 1(2) len by fastforce
              hence "⋀i. i < n ⟹ ratv2 $ i = gs.M.sumlist
                (map (λj. (of_int (B $$ (k, j)) ⋅⇩v (RAT fs) ! j)) [0..<k]) $ i"
                using fs.f_carrier 1(2) len v2carr gs.sumlist_nth sumlist_nth 
                  ratv2_def v2_def by simp
              then show ?thesis using ratv2carr carr by auto
            qed
            have "⋀i. i < k ⟹ (RAT fs) ! i = 
              gs.M.sumlist (map (λ j. gsi.μ i j ⋅⇩v gsi.gso j) [0 ..< Suc i])"
              using gsi.fi_is_sum_of_mu_gso len 1(2) by auto
            moreover have "⋀i. i < k ⟹ (λ j. gsi.μ i j ⋅⇩v gsi.gso j) ` {0 ..< Suc i}
              ⊆ gs.span (gsi.gso ` {0 ..<k})"
              using gs.span_mem gso0kcarr by auto
            ultimately have "⋀i. i < k ⟹ (RAT fs) ! i ∈ gs.span (gsi.gso ` {0 ..<k})"
              using gso0kcarr set_map set_upt vec_space.sumlist_in_span subsetD by smt
            then show ?thesis using rat atLeastLessThan_iff set_upt gso0kcarr imageE 
              set_map gs.smult_in_span vec_space.sumlist_in_span by smt
          qed
          have fs_gs_diff_span:
            "(RAT fs) !  k - fs'.gs.gso k ∈ gs.span (gsi.gso ` {0 ..<k})"
          proof -
            from fs'_gs_diff_span obtain v3 where sp: "v3 ∈ gs.span (gsi.gso ` {0 ..<k})"
              and eq: "(RAT fs) ! k - fs'.gs.gso k = v3 - ratv2" 
              using fs'.gs.gso_carrier len' 1(2) ratv2carr fs'id by fastforce
            have "v3+(-ratv2) ∈ gs.span(gsi.gso ` {0 ..<k})"
              by (metis sp gs.span_add1 gso0kcarr gram_schmidt.inv_in_span 
                  gso0kcarr ratv2span)
            moreover have "v3+(-ratv2) = v3-ratv2" using ratv2carr by auto
            ultimately have "v3 - ratv2 ∈ gs.span (gsi.gso ` {0 ..<k})" by simp
            then show ?thesis using eq by auto
          qed
          {
            fix i
            assume i: "i < k"
            hence "fs'.gs.gso k ∙ fs'.gs.gso i = 0" using 1(2) fs'.gs.orthogonal len' by auto
            hence "fs'.gs.gso k ∙ gsi.gso i = 0" using 1 i by simp
          }
          hence "⋀x. x ∈ gsi.gso ` {0..<k} ⟹ fs'.gs.gso k ∙ x = 0" by auto

          then show ?thesis
            using gsi.oc_projection_unique len len' fs_gs_diff_span 1(2) by auto
        qed
      qed
    qed

    have "⋀ i' j'. i' < m ⟹ j' < m ⟹ μ fs' i' j' = 
      (map_mat of_int (A i j c) * gsi.M m) $$ (i',j')" and
      "⋀ k. k < m ⟹ gso fs' k = gso fs k"
    proof -
      define rB where "rB = map_mat rat_of_int B"
      have rBcarr: "rB ∈ carrier_mat m m" using A(3) unfolding rB_def B_def by simp
      define rfs where "rfs = mat_of_rows n (RAT fs)"
      have rfscarr: "rfs ∈ carrier_mat m n" using Mfs unfolding rfs_def by simp

      {
        fix i'
        fix j'
        assume i': "i' < m"
        assume j': "j' < m"
        have prep: 
          "of_int_hom.vec_hom (row (B * mat_of_rows n fs) i') = row (rB * rfs) i'" 
          using len i' B_def A(3) rB_def rfs_def by (auto simp: scalar_prod_def)
        have prep2: "row (rB * rfs) i' = vec n (λl. row rB i' ∙ col rfs l)"
          using len fs.f_carrier i' B_def A(3) scalar_prod_def rB_def
          unfolding rfs_def by auto
        have prep3: "(vec m (λ j1. row rfs j1 ∙ gsi.gso j' / ∥gsi.gso j'∥⇧2)) =
          (vec m (λ j1. (gsi.M m) $$ (j1, j')))"
        proof -
          {
            fix x y
            assume x: "x < m" and y: "y < m"
            have "(gsi.M m) $$ (x,y) = (if y < x then map of_int_hom.vec_hom fs ! x 
              ∙ fs'.gs.gso y / ∥fs'.gs.gso y∥⇧2 else if x = y then 1 else 0)" 
              using gsi.μ.simps x y j' len gs_eq gsi.M_index by auto
            hence "row rfs x ∙ gsi.gso y / ∥gsi.gso y∥⇧2 = (gsi.M m) $$ (x,y)"
              unfolding rfs_def 
              by (metis carrier_matD(1) divide_eq_eq fs'.gs.β_zero fs'.gs.gso_norm_beta 
                  gs_eq gsi.μ.simps gsi.fi_scalar_prod_gso gsi.fs_carrier len len' 
                  length_map nth_rows rfs_def rfscarr rows_mat_of_rows x y)
          }
          then show ?thesis using j' by auto
        qed
        have prep4: "(1 / ∥gsi.gso j'∥⇧2) ⋅⇩v (vec m (λj1. row rfs j1 ∙ gsi.gso j')) =
          (vec m (λj1. row rfs j1 ∙ gsi.gso j' / ∥gsi.gso j'∥⇧2))" by auto

        have "map of_int_hom.vec_hom fs' ! i' ∙ fs'.gs.gso j' / ∥fs'.gs.gso j'∥⇧2
           = map of_int_hom.vec_hom fs' ! i' ∙ gsi.gso j' / ∥gsi.gso j'∥⇧2"
          using gs_eq j' by simp
        also have "… = row (rB * rfs) i' ∙ gsi.gso j' / ∥gsi.gso j'∥⇧2"
          using prep i' len' unfolding rB_def B_def by (simp add: fs'_prod)
        also have "… = 
          (vec n (λl. row rB i' ∙ col rfs l)) ∙ gsi.gso j' / ∥gsi.gso j'∥⇧2"
          using prep2 by auto
        also have "vec n (λl. row rB i' ∙ col rfs l) = 
            (vec n (λl. (∑j1=0..<m. (row rB i') $ j1 * (col rfs l) $ j1)))"
          using gsi.gso_carrier
          by (metis (no_types) carrier_matD(1) col_def dim_vec rfscarr scalar_prod_def)
        also have "… = 
            (vec n (λl. (∑j1=0..<m. rB $$ (i', j1) * rfs $$ (j1, l))))" 
          using rBcarr rfscarr i' by auto
        also have "… ∙ gsi.gso j' = 
            (∑j2=0..<n. (vec n 
            (λl. (∑j1=0..<m. rB $$ (i', j1) * rfs $$ (j1, l)))) $ j2 * (gsi.gso j') $ j2)"
          using gsi.gso_carrier len j' scalar_prod_def 
          by (smt gs.R.finsum_cong' gsi.gso_dim length_map)
        also have "… = (∑j2=0..<n.
            (∑j1=0..<m. rB $$ (i', j1) * rfs $$ (j1, j2)) * (gsi.gso j') $ j2)"
          using gsi.gso_carrier len j' by simp
        also have "… = (∑j2=0..<n. (∑j1=0..<m.
            rB $$ (i', j1) * rfs $$ (j1, j2) * (gsi.gso j') $ j2))" 
          by (smt gs.R.finsum_cong' sum_distrib_right)
        also have "… = (∑j1=0..<m. (∑j2=0..<n.
            rB $$ (i', j1) * rfs $$ (j1, j2) * (gsi.gso j') $ j2))"
          using sum.swap by auto
        also have "… = (∑j1=0..<m. rB $$ (i', j1) * (∑j2=0..<n. 
            rfs $$ (j1, j2) * (gsi.gso j') $ j2))"
          using gs.R.finsum_cong' sum_distrib_left by (smt gs.m_assoc)
        also have "… = row rB i' ∙ (vec m (λ j1. (∑j2=0..<n.
            rfs $$ (j1, j2) * (gsi.gso j') $ j2)))" 
          using rBcarr rfscarr i' scalar_prod_def
          by (smt atLeastLessThan_iff carrier_matD(1) carrier_matD(2) dim_vec 
              gs.R.finsum_cong' index_row(1) index_vec)
        also have "(vec m (λ j1. (∑j2=0..<n. rfs $$ (j1, j2) * (gsi.gso j') $ j2)))
            =  (vec m (λ j1. row rfs j1 ∙ gsi.gso j'))"
          using rfscarr gsi.gso_carrier len j' rfscarr by (auto simp add: scalar_prod_def)
        also have "row rB i' ∙ … / ∥gsi.gso j'∥⇧2 =
          row rB i' ∙ vec m (λ j1. row rfs j1 ∙ gsi.gso j' / ∥gsi.gso j'∥⇧2)"
          using prep4 scalar_prod_smult_right rBcarr carrier_matD(2) dim_vec row_def 
          by (smt gs.l_one times_divide_eq_left)
        also have "… = (rB * (gsi.M m)) $$ (i', j')" 
          using rBcarr i' j' prep3 gsi.M_def by (simp add: col_def)
        finally have 
          "map of_int_hom.vec_hom fs' ! i' ∙ fs'.gs.gso j' / ∥fs'.gs.gso j'∥⇧2 =
          (rB * (gsi.M m)) $$ (i', j')" by auto
      }
      then show "⋀ i' j'. i' < m ⟹ j' < m ⟹ μ fs' i' j' = 
        (map_mat of_int (A i j c) * gsi.M m) $$ (i',j')" 
        using B_def fs'.gs.β_zero fs'.gs.fi_scalar_prod_gso fs'.gs.gso_norm_beta
          len' rB_def by auto
      show "⋀ k. k < m ⟹ gso fs' k = gso fs k" using gs_eq by auto
    qed
  } note mu_gso = this

  show "⋀ k. k < m ⟹ gso fs' k = gso fs k" by fact
  {
    fix k
    have "k ≤ m ⟹ rat_of_int (d fs' k) = rat_of_int (d fs k)" for k
    proof (induct k)
      case 0
      show ?case by (simp add: d_def)
    next
      case (Suc k)
      hence k: "k ≤ m" "k < m" by auto 
      show ?case
        by (subst (1 2) LLL_d_Suc[OF _ k(2)], auto simp: Suc(1)[OF k(1)] mu_gso(2)[OF k(2)]
          LLL_invariant_weak_def lin lin' len len' latt latt')
    qed
    thus "k ≤ m ⟹ d fs' k = d fs k" by simp
  } note d = this
  {
    assume i': "i' < m" and j': "j' < m"
    have "μ fs' i' j' = (of_int_hom.mat_hom (A i j c) * gsi.M m) $$ (i',j')" by (rule mu_gso(1)[OF i' j'])
    also have "… = (if (i',j') = (i,j) then of_int c * gsi.d j else 0) + gsi.M m $$ (i',j')" 
      unfolding A(1) using i' j' by (auto simp: gsi.M_def)
    also have "gsi.M m $$ (i',j') = μ fs i' j'" 
      unfolding gsi.M_def using i' j' by simp
    also have "gsi.d j = of_int (d fs j)" 
      unfolding d_def by (subst Gramian_determinant_of_int[OF fs], insert ji i len, auto)
    finally show mu: "μ fs' i' j' = (if (i',j') = (i,j) then rat_of_int (c * d fs j) + μ fs i' j' else μ fs i' j')" 
      by simp
    let ?d = "d fs (Suc j')" 
    have d_fs: "of_int (dμ fs i' j') = rat_of_int ?d * μ fs i' j'" 
      unfolding dμ_def 
      using fs.fs_int_mu_d_Z_m_m[unfolded len, OF i' j'] 
      by (metis LLL.LLL.d_def assms(2) fs.fs_int_mu_d_Z_m_m fs_int.d_def i' 
          int_of_rat(2) j')
    have "rat_of_int (dμ fs' i' j') = rat_of_int (d fs' (Suc j')) * μ fs' i' j'" 
      unfolding dμ_def 
      using fs'.fs_int_mu_d_Z_m_m[unfolded len', OF i' j']
      using LLL.LLL.d_def fs'(1) fs'.dμ fs'.dμ_def fs_int.d_def i' j' by auto
    also have "d fs' (Suc j') = ?d" by (rule d, insert j', auto)
    also have "rat_of_int … * μ fs' i' j' = 
      (if (i',j') = (i,j) then rat_of_int (c * d fs j * ?d) else 0) + of_int (dμ fs i' j')" 
      unfolding mu d_fs by (simp add: field_simps)
    also have "… = rat_of_int ((if (i',j') = (i,j) then c * d fs j * ?d else 0) + dμ fs i' j')"
      by simp
    also have "… = rat_of_int ((if (i',j') = (i,j) then c * d fs j * d fs (Suc j) + dμ fs i' j' else dμ fs i' j'))"
      by simp
    finally show "dμ fs' i' j' = (if (i',j') = (i,j) then c * d fs j * d fs (Suc j) + dμ fs i' j' else dμ fs i' j')" 
      by simp
  }
qed

text ‹Eventually: Lemma 13 of Storjohann's paper.›
lemma mod_single_element: assumes lin: "lin_indep fs" 
  and len: "length fs = m" 
  and i: "i < m" and ji: "j < i"  
  and latt: "lattice_of fs = L" 
  and pgtz: "p > 0"
shows "∃ fs'. lattice_of fs' = L ∧ 
  map (map_vec (λ x. x mod p)) fs' = map (map_vec (λ x. x mod p)) fs ∧
  map (map_vec (λ x. x symmod p)) fs' = map (map_vec (λ x. x symmod p)) fs ∧
  lin_indep fs' ∧
  length fs' = m ∧ 
  (∀ k < m. gso fs' k = gso fs k) ∧ 
  (∀ k ≤ m. d fs' k = d fs k) ∧
  (∀ i' < m. ∀ j' < m. dμ fs' i' j' = (if (i',j') = (i,j) then dμ fs i j' symmod (p * d fs j' * d fs (Suc j')) else dμ fs i' j'))" 
proof -
  have inv: "LLL_invariant_weak fs" using LLL_invariant_weak_def assms by simp
  let ?mult = "d fs j * d fs (Suc j)" 
  define M where "M = ?mult" 
  define pM where "pM = p * M" 
  then have pMgtz: "pM > 0" using pgtz unfolding pM_def M_def using LLL_d_pos[OF inv] i ji by simp
  let ?d = "dμ fs i j" 
  define c where "c = - (?d symdiv pM)" 
  have d_mod: "?d symmod pM = c * pM + ?d" unfolding c_def using pMgtz sym_mod_sym_div by simp
  define A where "A = gram_schmidt_fs_int.inv_mu_ij_mat n (RAT fs)" 
  define fs' where fs': "fs' = Matrix.rows (A i j (c * p) * mat_of_rows n fs)"
  note main = change_single_element[OF lin len i ji A_def fs' latt]
  have "map (map_vec (λx. x mod p)) fs' = map (map_vec (λx. x mod p)) fs" 
    by (intro main, auto)
  from arg_cong[OF this, of "map (map_vec (poly_mod.inv_M p))"]
  have id: "map (map_vec (λx. x symmod p)) fs' = map (map_vec (λx. x symmod p)) fs" 
    unfolding map_map o_def sym_mod_def map_vec_map_vec .
  show ?thesis
  proof (intro exI[of _ fs'] conjI main allI impI id)
    fix i' j'
    assume ij: "i' < m" "j' < m" 
    have "dμ fs' i' j' = (if (i', j') = (i, j) then (c * p) * M + ?d else dμ fs i' j')" 
      unfolding main(8)[OF ij] M_def by simp
    also have "(c * p) * M + ?d = ?d symmod pM" 
      unfolding d_mod by (simp add: pM_def)
    finally show "dμ fs' i' j' = (if (i',j') = (i,j) then dμ fs i j' symmod (p * d fs j' * d fs (Suc j')) else dμ fs i' j')" 
      by (auto simp: pM_def M_def ac_simps)
  qed auto
qed 

text ‹A slight generalization to perform modulo on arbitrary set of indices $I$.›
lemma mod_finite_set: assumes lin: "lin_indep fs" 
  and len: "length fs = m" 
  and I: "I ⊆ {(i,j). i < m ∧ j < i}"
  and latt: "lattice_of fs = L" 
  and pgtz: "p > 0"
shows "∃ fs'. lattice_of fs' = L ∧
  map (map_vec (λ x. x mod p)) fs' = map (map_vec (λ x. x mod p)) fs ∧
  map (map_vec (λ x. x symmod p)) fs' = map (map_vec (λ x. x symmod p)) fs ∧
  lin_indep fs' ∧
  length fs' = m ∧ 
  (∀ k < m. gso fs' k = gso fs k) ∧ 
  (∀ k ≤ m. d fs' k = d fs k) ∧
  (∀ i' < m. ∀ j' < m. dμ fs' i' j' = 
     (if (i',j') ∈ I then dμ fs i' j' symmod (p * d fs j' * d fs (Suc j')) else dμ fs i' j'))"
proof -
  let ?exp = "λ fs' I i' j'. 
    dμ fs' i' j' = (if (i',j') ∈ I then dμ fs i' j' symmod (p * d fs j' * d fs (Suc j')) else dμ fs i' j')" 
  let ?prop = "λ fs fs'. lattice_of fs' = L ∧ 
    map (map_vec (λ x. x mod p)) fs' = map (map_vec (λ x. x mod p)) fs ∧
    map (map_vec (λ x. x symmod p)) fs' = map (map_vec (λ x. x symmod p)) fs ∧
    lin_indep fs' ∧
    length fs' = m ∧ 
    (∀ k < m. gso fs' k = gso fs k) ∧ 
    (∀ k ≤ m. d fs' k = d fs k)" 
  have "finite I" 
  proof (rule finite_subset[OF I], rule finite_subset)
    show "{(i, j). i < m ∧ j < i} ⊆ {0..m} × {0..m}" by auto
  qed auto
  from this I have "∃ fs'. ?prop fs fs' ∧ (∀ i' < m. ∀ j' < m. ?exp fs' I i' j')"
  proof (induct I)
    case empty
    show ?case
      by (intro exI[of _ fs], insert assms, auto)
  next
    case (insert ij I)
    obtain i j where ij: "ij = (i,j)" by force
    from ij insert(4) have i: "i < m" "j < i" by auto
    from insert(3,4) obtain gs where gs: "?prop fs gs" 
      and exp: "⋀ i' j'. i' < m ⟹ j' < m ⟹ ?exp gs I i' j'" by auto
    from gs have "lin_indep gs" "lattice_of gs = L" "length gs = m" by auto
    from mod_single_element[OF this(1,3) i this(2), of p] 
    obtain hs where hs: "?prop gs hs" 
      and exp': "⋀ i' j'. i' < m ⟹ j' < m ⟹ 
      dμ hs i' j' = (if (i', j') = (i, j) 
         then dμ gs i j' symmod (p * d gs j' * d gs (Suc j')) else dμ gs i' j')" 
      using pgtz by auto
    from gs i have id: "d gs j = d fs j" "d gs (Suc j) = d fs (Suc j)" by auto
    show ?case 
    proof (intro exI[of _ hs], rule conjI; (intro allI impI)?)
      show "?prop fs hs" using gs hs by auto
      fix i' j'
      assume *: "i' < m" "j' < m" 
      show "?exp hs (insert ij I) i' j'" unfolding exp'[OF *] ij using exp * i
        by (auto simp: id)
    qed
  qed
  thus ?thesis by auto
qed

end

end